GO Enrichment Analysis: Biological Processes
Francisca Samsing
18/12/2018
https://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/
# Load library and search for the organism
library(clusterProfiler)
library(pathview)
library(tidyverse)
library(enrichplot)Build the Salmon Data Base and import the gene universe (background genes from Salmon genome)
library(Category)
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, c("salmo salar","orgdb"))## AnnotationHub with 1 record
## # snapshotDate(): 2018-04-30
## # names(): AH61820
## # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
## # $species: Salmo salar
## # $rdataclass: OrgDb
## # $rdatadateadded: 2018-04-20
## # $title: org.Salmo_salar.eg.sqlite
## # $description: NCBI gene ID based annotations about Salmo salar
## # $taxonomyid: 8030
## # $genome: NCBI genomes
## # $sourcetype: NCBI/UniProt
## # $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.uniprot....
## # $sourcesize: NA
## # $tags: c("NCBI", "Gene", "Annotation")
## # retrieve record with 'object[["AH61820"]]'salmodb <- hub[["AH61820"]]
DatPkgFactory(salmodb)## An object of class "Org.XX.egDatPkg"
## Slot "name":
## [1] "OrgDb for Salmo salar"
##
## Slot "db":
## OrgDb object:
## | DBSCHEMAVERSION: 2.1
## | DBSCHEMA: NOSCHEMA_DB
## | ORGANISM: Salmo salar
## | SPECIES: Salmo salar
## | CENTRALID: GID
## | Taxonomy ID: 8030
## | Db type: OrgDb
## | Supporting package: AnnotationDbi
##
## Slot "installed":
## [1] FALSEcolumns(salmodb)## [1] "ACCNUM" "ALIAS" "CHR" "ENTREZID" "EVIDENCE"
## [6] "EVIDENCEALL" "GENENAME" "GID" "GO" "GOALL"
## [11] "ONTOLOGY" "ONTOLOGYALL" "PMID" "REFSEQ" "SYMBOL"
## [16] "UNIGENE"gene_universe <- read.csv('Results/gene_universe.csv')PATHWAY ANALYSIS FOR POMV 6 HPI
The following code loads ALL differentially expressed genes with an adjusted p-value < 0.05, and then selects those with an absolute logfold change of 2.
POMV6 <- read.csv('Results/ControlvsPOMV6_ALL.csv') # CHECK FOR NAs
POMV6_SigGeneList <- POMV6[abs(POMV6$logFC) >= 2,]
rownames(POMV6_SigGeneList) <- NULL
POMV6_SigGeneList_logFC <- POMV6_SigGeneList$logFC
names(POMV6_SigGeneList_logFC) <- POMV6_SigGeneList$ENTREZID
ego_POMV6 <- enrichGO(gene = names(POMV6_SigGeneList_logFC),
universe = as.character(gene_universe$ENTREZID),
keyType = 'ENTREZID',
OrgDb = salmodb,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
head(summary(ego_POMV6))## ID Description GeneRatio
## GO:0006952 GO:0006952 defense response 14/69
## GO:0002376 GO:0002376 immune system process 17/69
## GO:0009615 GO:0009615 response to virus 7/69
## GO:0006955 GO:0006955 immune response 13/69
## GO:0051607 GO:0051607 defense response to virus 6/69
## GO:0043207 GO:0043207 response to external biotic stimulus 7/69
## BgRatio pvalue p.adjust qvalue
## GO:0006952 132/10679 8.378159e-14 2.211834e-11 2.010758e-11
## GO:0002376 277/10679 1.186589e-12 1.566297e-10 1.423906e-10
## GO:0009615 29/10679 4.799528e-10 4.223585e-08 3.839623e-08
## GO:0006955 221/10679 1.247334e-09 8.232402e-08 7.484002e-08
## GO:0051607 21/10679 2.930697e-09 1.547408e-07 1.406734e-07
## GO:0043207 42/10679 7.764224e-09 2.928222e-07 2.662020e-07
## geneID
## GO:0006952 100195148/100196303/100136541/100136587/106566099/106570878/106570889/106591927/106596334/106600865/100194889/100136920/100196194/106602560
## GO:0002376 100195148/100196303/100136541/100136587/106565339/106565855/106566099/106570878/106570889/106591927/106596334/106600865/100194889/100136920/100196194/106602560/106586517
## GO:0009615 100196303/100136541/106566099/106596334/106600865/100194889/106602560
## GO:0006955 100195148/100196303/100136541/100136587/106565339/106565855/106570878/106570889/106591927/100136920/100196194/106602560/106586517
## GO:0051607 100136541/106566099/106596334/106600865/100194889/106602560
## GO:0043207 100196303/100136541/106566099/106596334/106600865/100194889/106602560
## Count
## GO:0006952 14
## GO:0002376 17
## GO:0009615 7
## GO:0006955 13
## GO:0051607 6
## GO:0043207 7barplot(ego_POMV6, showCategory=20, font.size = 10)
dotplot(ego_POMV6, showCategory=20, font.size = 10)
emapplot(ego_POMV6)
cnetplot(ego_POMV6, categorySize="genNUM", foldChange=POMV6_SigGeneList_logFC)
PATHWAY ANALYSIS FOR POMV 24 HPI
The following code loads ALL differentially expressed genes with an adjusted p-value < 0.05, and then selects those with an absolute logfold change of 2.
POMV24 <- read.csv('Results/ControlvsPOMV24_ALL.csv') # CHECK FOR NAs
POMV24_SigGeneList <- POMV24[abs(POMV24$logFC) >= 2,]
rownames(POMV24_SigGeneList) <- NULL
POMV24_SigGeneList_logFC <- POMV24_SigGeneList$logFC
names(POMV24_SigGeneList_logFC) <- POMV24_SigGeneList$ENTREZID
ego_POMV24 <- enrichGO(gene = names(POMV24_SigGeneList_logFC),
universe = as.character(gene_universe$ENTREZID),
keyType = 'ENTREZID',
OrgDb = salmodb,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
head(summary(ego_POMV24))## ID Description GeneRatio
## GO:0006952 GO:0006952 defense response 55/1356
## GO:0002376 GO:0002376 immune system process 84/1356
## GO:0006955 GO:0006955 immune response 70/1356
## GO:0051607 GO:0051607 defense response to virus 17/1356
## GO:0007264 GO:0007264 small GTPase mediated signal transduction 96/1356
## GO:0098542 GO:0098542 defense response to other organism 20/1356
## BgRatio pvalue p.adjust qvalue
## GO:0006952 132/10679 5.941100e-17 6.951087e-14 6.666539e-14
## GO:0002376 277/10679 3.847759e-15 2.250939e-12 2.158795e-12
## GO:0006955 221/10679 7.496142e-14 2.923496e-11 2.803820e-11
## GO:0051607 21/10679 1.918880e-12 5.612724e-10 5.382964e-10
## GO:0007264 394/10679 8.178118e-11 1.913680e-08 1.835342e-08
## GO:0098542 33/10679 1.155639e-10 2.253497e-08 2.161248e-08
## geneID
## GO:0006952 100196184/106601862/100195999/100195148/100137019/100136436/100195681/100136439/100136449/100136453/100136456/100136541/100136587/100302030/101448043/106561958/106565979/106566099/106569524/106570878/106570889/106574371/106575194/106576643/106576647/106577834/106577950/106583498/106589772/106590949/106591927/106593218/106595807/106596088/106596190/106596334/106599427/106600865/106600963/106600969/106600975/106601371/106601388/106604275/106606549/106607463/106609132/106609858/106611470/100194889/100136920/100196194/100270812/106602560/106611198
## GO:0002376 100196184/106601862/100195999/100195148/100329176/100137019/100136436/100195681/100136548/100846970/106584118/100136449/100136541/100136587/100302030/101448043/106563358/106564720/106565339/106565694/106565799/106565855/106565979/106566099/106567717/106569524/106570878/106570886/106570889/106572809/106574371/106574975/106575194/106575840/106575849/106576643/106576647/106577010/106577763/106577765/106577834/106577950/106579494/106581031/106585878/106586893/106588814/106589772/106590519/106590949/106591927/106593218/106595807/106596088/106596334/106599427/106600865/106600963/106600975/106601371/106601388/106602980/106604275/106605663/106606549/106606848/106607463/106607531/106608811/106609048/106609132/106609531/106609858/106611470/106613961/100194889/100136920/100196194/100270812/106602560/100136458/106586517/100196575/106611198
## GO:0006955 100196184/106601862/100195999/100195148/100137019/100136436/100195681/100136548/100846970/106584118/100136449/100136541/100136587/100302030/106563358/106564720/106565339/106565694/106565799/106565855/106565979/106567717/106569524/106570878/106570886/106570889/106572809/106574371/106574975/106575194/106575840/106575849/106576643/106576647/106577010/106577765/106577834/106577950/106579494/106581031/106585878/106586893/106589772/106590519/106591927/106593218/106595807/106596088/106599427/106600975/106601371/106601388/106602980/106604275/106605663/106606549/106606848/106607531/106608811/106609048/106609858/106613961/100136920/100196194/100270812/106602560/100136458/106586517/100196575/106611198
## GO:0051607 100137019/100136436/100136541/100302030/101448043/106566099/106590949/106596334/106600865/106600963/106600975/106607463/106609132/106611470/100194889/106602560/106611198
## GO:0007264 106574446/106610330/100380466/106613922/100196534/106560488/106560490/106561157/106562829/106562943/106563223/106565273/106565544/106565828/106566806/106566812/106566904/106567797/106568132/106568270/106568392/106569172/106570905/106571247/106571902/106572387/106576073/106576227/106576952/106577639/106577853/106577869/106578971/106579157/106579284/106579856/106579924/106579959/106580138/106580155/106580158/106580556/106581424/106581831/106582025/106583366/106584026/106584283/106584386/106584876/106585456/106585561/106585588/106586158/106586233/106588489/106588797/106588976/106589836/106591236/106591860/106591994/106592352/106594155/106594609/106599451/106600087/106600822/106601059/106601313/106603123/106603740/106604741/106605732/106606050/106607078/106607387/106607499/106608153/106608271/106608981/106609016/106609220/106609966/106611722/106611922/106612426/106613940/106581351/106601189/106610804/106613823/106564550/106600638/106611626/106607670
## GO:0098542 106601862/100137019/100136436/100136439/100136453/100136541/100302030/101448043/106566099/106590949/106596334/106600865/106600963/106600975/106607463/106609132/106611470/100194889/106602560/106611198
## Count
## GO:0006952 55
## GO:0002376 84
## GO:0006955 70
## GO:0051607 17
## GO:0007264 96
## GO:0098542 20barplot(ego_POMV24, showCategory=20, font.size = 10)
dotplot(ego_POMV24, showCategory=20, font.size = 10)
emapplot(ego_POMV24, showCategory = 20)
cnetplot(ego_POMV24, categorySize="genNUM", foldChange=POMV24_SigGeneList_logFC, node_label = FALSE)
PATHWAY ANALYSIS FOR ISAV 6 HPI
The following code loads ALL differentially expressed genes with an adjusted p-value < 0.05, and then selects those with an absolute logfold change of 2.
ISAV6 <- read.csv('Results/ControlvsISAV6_ALL.csv') # CHECK FOR NAs
ISAV6_SigGeneList <- ISAV6[abs(ISAV6$logFC) >= 2,]
rownames(ISAV6_SigGeneList) <- NULL
ISAV6_SigGeneList_logFC <- ISAV6_SigGeneList$logFC
names(ISAV6_SigGeneList_logFC) <- ISAV6_SigGeneList$ENTREZID
ego_ISAV6 <- enrichGO(gene = names(ISAV6_SigGeneList_logFC),
universe = as.character(gene_universe$ENTREZID),
keyType = 'ENTREZID',
OrgDb = salmodb,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
head(summary(ego_ISAV6))## ID Description GeneRatio BgRatio
## GO:0002376 GO:0002376 immune system process 27/85 277/10679
## GO:0006955 GO:0006955 immune response 24/85 221/10679
## GO:0006952 GO:0006952 defense response 18/85 132/10679
## GO:0045087 GO:0045087 innate immune response 11/85 71/10679
## GO:0034097 GO:0034097 response to cytokine 8/85 37/10679
## GO:0050776 GO:0050776 regulation of immune response 9/85 70/10679
## pvalue p.adjust qvalue
## GO:0002376 1.219522e-22 3.243928e-20 2.323510e-20
## GO:0006955 3.202610e-21 4.259471e-19 3.050907e-19
## GO:0006952 8.325678e-18 7.382101e-16 5.287536e-16
## GO:0045087 7.253815e-12 4.823787e-10 3.455106e-10
## GO:0034097 3.684795e-10 1.960311e-08 1.404101e-08
## GO:0050776 3.645326e-09 1.616095e-07 1.157551e-07
## geneID
## GO:0002376 100195999/100195148/100196303/100136548/100846970/100136541/100136587/106563358/106564720/106565339/106565855/106566099/106570878/106570889/106574975/106575849/106585878/106591927/106593218/106596088/106599427/106600865/100194889/100136920/100196194/100270812/106602560
## GO:0006955 100195999/100195148/100196303/100136548/100846970/100136541/100136587/106563358/106564720/106565339/106565855/106570878/106570889/106574975/106575849/106585878/106591927/106593218/106596088/106599427/100136920/100196194/100270812/106602560
## GO:0006952 100195999/100195148/100196303/100136541/100136587/106566099/106570878/106570889/106591927/106593218/106596088/106599427/106600865/100194889/100136920/100196194/100270812/106602560
## GO:0045087 100195148/100196303/100136541/100136587/106591927/106593218/106596088/106599427/100136920/100270812/106602560
## GO:0034097 100136541/100136587/106591927/106593218/106596088/106599427/100136920/100270812
## GO:0050776 100136548/100846970/106574975/106575849/106591927/106593218/106596088/106599427/106602560
## Count
## GO:0002376 27
## GO:0006955 24
## GO:0006952 18
## GO:0045087 11
## GO:0034097 8
## GO:0050776 9barplot(ego_ISAV6, showCategory=20, font.size = 10)
dotplot(ego_ISAV6, showCategory=20, font.size = 10)
emapplot(ego_ISAV6, font.size = 10)
cnetplot(ego_ISAV6, categorySize="genNUM", foldChange=POMV6_SigGeneList_logFC, font.size = 10)
PATHWAY ANALYSIS FOR ISAV 24 HPI
The following code loads ALL differentially expressed genes with an adjusted p-value < 0.05, and then selects those with an absolute logfold change of 2.
ISAV24 <- read.csv('Results/ControlvsISAV24_ALL.csv') # CHECK FOR NAs
ISAV24_SigGeneList <- ISAV24[abs(ISAV24$logFC) >= 2,]
rownames(ISAV24_SigGeneList) <- NULL
ISAV24_SigGeneList_logFC <- ISAV24_SigGeneList$logFC
names(ISAV24_SigGeneList_logFC) <- ISAV24_SigGeneList$ENTREZID
ego_ISAV24 <- enrichGO(gene = names(ISAV24_SigGeneList_logFC),
universe = as.character(gene_universe$ENTREZID),
keyType = 'ENTREZID',
OrgDb = salmodb,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
head(summary(ego_ISAV6))## ID Description GeneRatio BgRatio
## GO:0002376 GO:0002376 immune system process 27/85 277/10679
## GO:0006955 GO:0006955 immune response 24/85 221/10679
## GO:0006952 GO:0006952 defense response 18/85 132/10679
## GO:0045087 GO:0045087 innate immune response 11/85 71/10679
## GO:0034097 GO:0034097 response to cytokine 8/85 37/10679
## GO:0050776 GO:0050776 regulation of immune response 9/85 70/10679
## pvalue p.adjust qvalue
## GO:0002376 1.219522e-22 3.243928e-20 2.323510e-20
## GO:0006955 3.202610e-21 4.259471e-19 3.050907e-19
## GO:0006952 8.325678e-18 7.382101e-16 5.287536e-16
## GO:0045087 7.253815e-12 4.823787e-10 3.455106e-10
## GO:0034097 3.684795e-10 1.960311e-08 1.404101e-08
## GO:0050776 3.645326e-09 1.616095e-07 1.157551e-07
## geneID
## GO:0002376 100195999/100195148/100196303/100136548/100846970/100136541/100136587/106563358/106564720/106565339/106565855/106566099/106570878/106570889/106574975/106575849/106585878/106591927/106593218/106596088/106599427/106600865/100194889/100136920/100196194/100270812/106602560
## GO:0006955 100195999/100195148/100196303/100136548/100846970/100136541/100136587/106563358/106564720/106565339/106565855/106570878/106570889/106574975/106575849/106585878/106591927/106593218/106596088/106599427/100136920/100196194/100270812/106602560
## GO:0006952 100195999/100195148/100196303/100136541/100136587/106566099/106570878/106570889/106591927/106593218/106596088/106599427/106600865/100194889/100136920/100196194/100270812/106602560
## GO:0045087 100195148/100196303/100136541/100136587/106591927/106593218/106596088/106599427/100136920/100270812/106602560
## GO:0034097 100136541/100136587/106591927/106593218/106596088/106599427/100136920/100270812
## GO:0050776 100136548/100846970/106574975/106575849/106591927/106593218/106596088/106599427/106602560
## Count
## GO:0002376 27
## GO:0006955 24
## GO:0006952 18
## GO:0045087 11
## GO:0034097 8
## GO:0050776 9barplot(ego_ISAV24, showCategory=20, font.size = 10)
dotplot(ego_ISAV24, showCategory=20, font.size = 10)
emapplot(ego_ISAV24, font.size = 10)
cnetplot(ego_ISAV24, categorySize="genNUM", foldChange=POMV6_SigGeneList_logFC, font.size = 10)